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Abstract: Lattice simulations have demonstrated that a background (electro)magnetic field reduces 
the chiral/deconfinement transition temperature of quantum chromodynamics for eB < 1 GeV^. On 
the level of observables, this reduction manifests itself in an enhancement of the Polyakov loop and 
in a suppression of the light quark condensates (inverse magnetic catalysis) in the transition region. 
In this paper, we report on lattice simulations of 1 + 1 + 1-flavor QCD at an unprecedentedly high 
value of the magnetic field eB = 3.25 GeV^. Based on the behavior of various observables, it is shown 
that even at this extremely strong field, inverse magnetic catalysis prevails and the transition, albeit 
becoming sharper, remains an analytic crossover. In addition, we develop an algorithm to directly 
simulate the asymptotically strong magnetic field limit of QCD. We find strong evidence for a first- 
order deconfinement phase transition in this limiting theory, implying the presence of a critical point 
in the QCD phase diagram. Based on the available lattice data, we estimate the location of the critical 
point. 
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1 Introduction 

Quantum chromodynamics (QCD) exhibits a finite temperature transition separating the chirally bro¬ 
ken, low-temperature phase from the chirally symmetric, high-temperature regime, where quarks and 
gluons are deconfined. Although this transition is no real phase transition but merely an analytic 
crossover [1], it is marked by the pronounced behavior of the corresponding (approximate) order pa¬ 
rameters: the drop in the light quark condensates, accompanied by the increase in the Polyakov loop. 
The characteristic dependence of the transition temperature on further parameters of the system probes 
our understanding of QCD and maps out the phase diagram in the corresponding parameter space. One 
parameter that is thought to have rich physical applications, ranging from neutron star physics through 
heavy-ion collisions to the cosmology of the early universe, is a background (electro)magnetic field 
B. The relevant range of magnetic fields, where QCD interactions compete with the electromagnetic 
forces, is given by multiples of the pion mass squared m^. We refer the reader to recent reviews on the 
subject [2-4]. 

QCD with background magnetic fields can be studied directly using non-perturbative lattice sim¬ 
ulations. Continuum extrapolated results employing staggered quarks with physical masses have been 
used to map out the phase diagram [5] for 0 < eB < 1 GeV^. According to these results, the magnetic 
field increases the light quark condensates well below and well above Tf. (magnetic catalysis) but de¬ 
creases them in the transition region [6] (inverse magnetic catalysis). As a result of this non-monotonous 
dependence of the condensate on B and on T, the transition temperature is significantly reduced by 
the magnetic field. The same tendency has been observed for the Polyakov loop as well, giving a sim¬ 
ilarly decreasing transition temperature [7]. More recent lattice simulations employing different quark 
discretizations are consistent with this picture [8]. 

The magnetic catalysis of the condensate at low temperatures is a very robust concept. It arises 
naturally due to the Landau-level structure of charged particle energies in the presence of B. In the 
strong field limit, magnetic catalysis can be understood in terms of the dimensional reduction of the 
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system and the high degeneracy of the lowest Landau-level [9, 10]. For low magnetic fields, it can 
be related to the positivity of the QED /3-function that fixes the dependence of the condensate on 
B to order [11, 12], Magnetic catalysis even has connections to solid state physics models like the 
Hofstadter model [13]. In line with these arguments, the catalysis of the condensate at low temperatures 
was observed in a variety of model settings and effective theories of QCD. However, in most of these 
models, magnetic catalysis takes place not only for T < Tc, but for all temperatures, giving rise to a 
monotonously increasing dependence of Tc on B. Thus, for the phase diagram, these models predict 
just the opposite of what the above discussed lattice results suggest. For a recent summary on these 
model approaches and a comparison to the lattice results, see Refs. [3, 14]. 

While the mechanisms behind magnetic catalysis, as mentioned above, are quite transparent, the 
opposite behavior around Tc - inverse magnetic catalysis - apparently has its origin in the rearrangement 
of the gluonic configurations that dominate the QCD path integral and is thus highly nontrivial [7]. Sev¬ 
eral attempts have been made recently to understand this behavior in effective approaches to QCD [15], 
among others, by introducing new, H-dependent model parameters or by taking into account the run¬ 
ning of the QCD coupling with the magnetic field. Several of these models exhibit a non-monotonous 
Tc{B) dependence, with an initial reduction followed by an enhancement due to the magnetic field. In 
certain settings, it was even shown that no matter how the existing parameters of the model are tuned 
as functions of 73, the transition temperature always tends to rise above a given threshold magnetic 
field [16]. 

It is just the apparent universality of magnetic catalysis that has made the lattice results about 
inverse catalysis and the decreasing Tc{B) dependence for 0 < eB < 1 GeV^ so unexpected. It was 
speculated that magnetic catalysis should reappear at even stronger magnetic fields, and different 
hypotheses were recently put forward about the strong B regime of the phase diagram. In particular, 
the strong B limit was argued to induce a new critical point [17]. The transition temperature was 
conjectured to turn around and increase if the magnetic field is sufficiently strong [18-21]. In other 
cases, Tc was argued to keep decreasing and to hit zero [22]. The transition temperatures for chiral 
restoration and for deconfinement were predicted to split [23] in the presence of the magnetic field, and 
a splitting between the chiral restoration temperature for the up and down quarks was also argued to 
take place [21, 24]. Let us refer the reader to the recent reviews [3, 14] for details. 

In this paper, we aim to check these conjectures by means of first-principles lattice simulations of 
1 -|- 1 -|- 1-flavor QCD at an unprecedentedly strong magnetic field eB = 3.25 GeV^. In addition, we 
also simulate the 73 —>■ oo limit directly, by considering the effective theory relevant for this limit [25]. 
We find strong evidence that this limiting theory has a first-order deconfinement phase transition and, 
thus, the QCD phase diagram exhibits a critical point at strong magnetic fields, where the analytic 
crossover terminates. Based on our results, we estimate the location of the critical point, and sketch the 
dependence of the deconfinement transition temperature on 73 over a broad range. Besides answering 
a fundamental question about the QCD phase diagram, we believe that the results will also be useful 
for building/refining effective theories and models of QCD. 

The paper is organized as follows. In Sec. 2 we discuss the details of the simulations and define 
the observables used to study the phase diagram. Sec. 3 contains the lattice results in ordinary QCD, 
followed by Sec. 4, where we discuss the simulations of the anisotropic theory in the asymptotic limit. 
The derivation of this effective theory and the employed simulation algorithms are discussed in the 
appendices App. A and App. B. Finally, in Sec. 5 we summarize our findings regarding the QCD phase 
diagram and conclude. 
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2 Setup and observables 


We consider a spatially symmetric iV| x Nt lattice with spacing a so that the temperature is given 
by T = the spatial volume hy V = (Nga)^ and the four-volume by V 4 = V/T. Given this 

geometry, we simulate 1 -|- 1 -|- 1-flavor QCD, described by the partition function, 

Z= [det (2.1) 

f=u,d,s 

given by the functional integral over the gluonic links U. We employ stout smeared rooted staggered 
quarks described by the fermion matrix M. In Eq. (2.1), Sg = SgV^ is the tree-level Symanzik improved 
gauge action and f3 = dj^ the inverse gauge coupling. For further details of the simulation setup and 
algorithm, see Refs. [5, 26]. The parameters of the fermion matrix are the quark masses = rrid rUg 
and the electric charges Qu = —2qd = —2qs = 2e/3 (e > 0 is the elementary charge), which enter 
in the product with B. The quark masses are set to their physical values along the line of constant 
physics [27]. The magnetic field is oriented along the positive z-direction and has the quantized flux 

= {aNsf eB = GirNb, G 0 < Nb < Ng, (2.2) 


where we used that the smallest charge in the system is that of the down quark. Lattice discretization 
effects are suppressed as long as the flux quantum Nb is much smaller than the period N^. Previous 
experience suggests that Nb < Ng /IG is a reasonable choice [5]. 

In the following, we employ the fixed-W approach and change the temperature T{/3) = (W a(/3))“^ 
by varying the inverse gauge coupling. This also implies that a given flux quantum corresponds to 
different magnetic fields at different temperatures, i.e. eB oc NbT‘^{j3). In particular, we choose (3 
values where a fixed magnetic field eB = 3.25 GeV^ is represented by integer flux quanta. Although 
this implies that only discrete temperatures are allowed, at this strong magnetic field the temperature 
differences are small enough in order to map out the transition region (see below). 

Next, we define the observables that can be used to pin down the transition temperature. We begin 
with the quark condensates and susceptibilities, signaling chiral symmetry. 




1 d\ogZ 
V 4 dmf 


iXf) 


druf ’ 


(2.3) 


and employ the normalization inspired by the Gell-Mann-Oakes-Renner relation, introduced in Ref. [6] 
for the condensate. 
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Here, = 135 MeV is the pion mass and F = 8 G MeV the chiral limit of the pion decay constant. 
Both Tj and ^-re free of additive and of multiplicative divergences [5]. In addition, S is normalized 
to be unity at T = R = 0 and approaches zero above the transition region [6]. The vacuum values 
necessary for the additive renormalization were determined in Refs. [5, 6]. 

The approximate order parameter for center symmetry, related to the deconfinement transition, is 
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the Polyakov loop, defined on the lattice as 


P=i^^TrnC/4(t,x)^. (2.5) 

In full QCD, the fermion determinant breaks center symmetry explicitly, so that the spontaneous 
breaking always occurs towards the real center element and (ReP) is a valid (approximate) order 
parameter. In pure gauge theory (this will be relevant for the P —)• oo limit, see Sec. 4), there is no 
explicit breaking and the three center sectors are equivalent. In this case, it is convenient to consider 
the projection of the Polyakov loop to the nearest center element (see, e.g.. Ref. [28]), 


ppr 


ReP, 

< Re[Pe-*2p3]^ 
Re[Pe*2vr/3j^ 


argP e [-7r/3,7r/3], 
argP e (7r/3,7r], 
argP G (—vr, —7r/3). 


( 2 . 6 ) 


Simulating pure gauge theory on a finite lattice, (P) always vanishes due to the tunneling between 
center sectors, while (PP'") is positive in the deconfined phase. The susceptibility of the projected 
Polyakov loop is defined as 


Xppr = V 


^ppr2) _ 


(2.7) 


The Polyakov loop renormalizes multiplicatively, with a temperature-dependent renormalization 
constant 

Pr{T,B) = Z{T)- P{T,B) (2.8) 


which is determined by enforcing (Pr(r*,0)) = P* and we chose T* = 162 MeV and P* = 1. The 
renormalization was discussed in detail and Z(T) was determined in Ref. [7]. Notice that while ReP < 3 
by construction, the renormalized observable has no upper bound. 

An observable that strongly correlates with P - and, thus, is sensitive to the deconfinement tran¬ 
sition - is the strange quark number susceptibility. 


1 (9^ log Z 

V4 ■ r2 5^2 


(2.9) 


Note that c| contains neither additive nor multiplicative divergences. 
Finally, the trace anomaly 

j(4>) ^ _ L 


V 4 d log a 


can be written as a sum of gluonic, fermionic and magnetic contributions [ 12 ], 


( 2 . 10 ) 


li^)(B,T) 


dl3 d{mfa) 

5 log a ^ (9 logo 


+ bi{eBf, 


( 2 . 11 ) 


where 61 = X]j( 9 //e)^/( 47 r^) is the lowest-order QED /3-function coefficient. The magnetic term ap¬ 
pears due to electric charge renormalization and stems from the counter-term canceling the P-dependent 
additive divergence of the thermodynamic potential logZ [12]. Note that this term is finite and in¬ 
dependent of the regularization, once the continuum limit is taken, see discussion in Ref. [12]. Notice 
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furthermore that the derivative in the dehnition of I is evaluated at hxed magnetic flux and not at 
hxed magnetic held eB [this is indicated by the superscript ($)]. The need for distinguishing between 
the two directional derivatives was hrst discussed in Ref. [29] and put into practice for the trace anomaly 
in Ref. [12]. 

3 Results at eB = 3.25 GeV^ 

We extend the previously published data on the light condensates and susceptibilities [6], on the 
Polyakov loop [7], on the strange quark number susceptibility [5], and on the trace anomaly [12, 29] 
using our new results at eB = 3.25 GeV^. To achieve this magnetic held strength, a temporal lattice 
extent Nt = 16 turned out to be necessary. These Nt = 16 lattices are hner than the hnite temperature 
conhgurations used in Refs. [5-7] (W = 6, 8 and 10) to extrapolate to the continuum limit. Thus, our 
results - although not strictly continuum extrapolated - are expected to lie close to the limit a —?• 0. 
We use two spatial lattice sizes 32^ x 16 and 48^ x 16 to control hnite size effects. 

Let us start the discussion with the light quark condensates. The average of and is plotted 
in the left panel of Fig. 1, compared to the B = 0 and B = 1 GeV^ continuum extrapolated results [6]. 
In addition to the data at nonzero temperatures, we also indicate an estimate for the zero-temperature 
condensate. This is obtained by htting and extrapolating the available lattice data at T = 0 by a 
free-theory inspired form ~ BlogB. The systematic uncertainty is taken into account by varying the 
ht interval. 

While the condensate is increased by the magnetic held at low temperatures, rehecting the well- 
known magnetic catalysis effect, the results also clearly show the reduction of in the transition 

region. Thus, inverse magnetic catalysis is observed to persist in the transition region even for our 
strong magnetic held eB = 3.25 GeV^, pushing Tc further down. In particular, we employ the inhection 
point of the average condensate to hnd Tc{'Sud} = 112(3) MeV. As a side-remark, we mention that 
since the transition region is shifted to considerably lower temperatures, the vacuum values determined 




T (MeV) 1 (l^eV) 


Figure 1. Left panel: average light quark condensate as a function of the temperature for three different 
magnetic ftelds. Right panel: up (open points) and down (ftlled points) quark condensates for the same set of 
magnetic ftelds. The curves are spline interpolations and merely serve to guide the eye. 
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Figure 2. Left panel: the Polyakov loop for three values of eB. At the highest magnetic field, the curve is a 
spline interpolation, while for the lower fields the band is the result of a combined continuum extrapolation and 
interpolation in T [7]. Right panel: the strange quark number susceptibility for the same set of magnetic fields. 
For the highest magnetic field, a spline interpolation is shown, whereas for the lower fields the bands represent 
a continuum estimate based on the results of Ref. [5]. 

for 3.45 < (5 < 3.85 in Refs. [5, 6] suffice to perform the additive renormalization of the condensates, 
and there is no need for additional T = 0 simulations on finer lattices. 

Due to the different electric charges, is expected to be more sensitive to the magnetic field 
than Tid- On that account, even a splitting in the transition temperatures might seem plausible, see 
Refs. [21, 24]. To check whether this is the case, in the right panel of Fig. 1 the two condensates are 
plotted separately. Even though the difference is pronounced throughout the temperature 

range in question, fitting for the inffection points gives consistent values Tc{S„} = 112(3) MeV and 
Tc{S^} = 111(3) MeV. An apparent implication of this finding is that the temperature, at which 
the transition between the chirally broken and restored phases takes place, is encoded in the gluonic 
configurations rather than in the operator insertion. In lattice language; Tc seems to be a quantity driven 
predominantly by sea and not by valence effects. This also suggests that purely gluonic observables 
would also exhibit similar transition temperatures. 

This brings us to the simplest, purely gluonic quantity: the Polyakov loop (2.5). The (real part of 
the) renormalized observable is plotted in the left panel of Fig. 2, for the same set of magnetic fields, and 
is observed to be drastically enhanced by the magnetic field for all temperatures. The inflection point 
of Pr{T) is much more pronounced as compared to the case at = 1 GeV^ and is determined to be 
Tc{P} = 109(3) MeV. This value is indeed consistent with the transition temperatures obtained above 
for the light quark condensates. We conclude that the gluonic configurations are vastly different on the 
two ‘sides’ of the transition, and predestine the behavior of the light condensates, independently of the 
electric charge that appears in the operator. We also observe the strange quark number susceptibility 
to exhibit an analogous trend, see the right panel of Fig. 2. Performing a similar fit as for Pj., we obtain 
Tc{c|} = 109(3) MeV. 

A further observable of interest for the QCD equation of state is the trace anomaly (2.10). It mea¬ 
sures the breaking of conformal symmetry by the gluonic condensate, by the quark condensates and by 
the magnetic field itself. As B grows, the latter effect becomes dominant and is increased drastically, 
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as visible in Fig. 3. Since contains in¬ 
dependent contributions already at zero tem¬ 
perature, the usual normalization fT‘^ is 
not useful [12], This large T = 0 contribution 
also damps the behavior of in the transi¬ 
tion region. The small kink around Tc, moving 
towards smaller temperatures as B grows, is 
to some extent still visible. We note that in 
order to determine further quantities related 
to the equation of state (e.g. pressure, entropy 
density etc.), one would need additional simu¬ 
lations at low temperature (see the method de¬ 
veloped in Ref. [12]). This is outside the scope 
of the present paper. 

Besides the characteristic temperature, the 
strength of the transition at high magnetic fields is also of interest. To determine, whether the smooth 
crossover at eB < 1 GeV^ turns into a real phase transition at eB = 3.25 GeV^, we analyze the 
average of the light quark susceptibilities Xu Xd- This observable exhibits a peak at the transition 
temperature, see the left panel of Fig. 4. For real phase transitions, the height h of this peak diverges in 
the infinite volume limit: h oc V for first-order transitions and h oc with a critical exponent a < 1 
for second-order transitions. In contrast, for the case of an analytic crossover, h is independent of the 
volume. To perform this finite size scaling study, we compare the results obtained on the 48^ x 16 and 
on the 32^ x 16 ensembles. The left panel of Fig. 4 shows no sign of a singularity as V is increased (note 
that for a first-order transition, the peak heights for the two volumes would differ by more than three). 
This leads us to conclude that the transition remains an analytic crossover even at eB = 3.25 GeV^. 
We mention moreover that finite volume effects are also absent from the other observables discussed 
above. 

Although the transition remains an analytic crossover, it is instructive to analyze the R-dependence 
of the susceptibility peak in more detail. We normalize Xu + such that its peak maximum equals 





Figure 3. The trace anomaly for three different magnetic 
fields. Note the logarithmic scale. 


Figure 4. Left panel: finite size scaling of the average light quark susceptibility. Right panel: the dependence 
of the peak width on the magnetic field. 

















unity, and plot it in the right panel of Fig. 4 against the temperature. Here, T is shifted so that the 
observable equals 0.5 at zero. Then, the peak width w{B) at half maximum can be read off at the 
rightmost intersection of the observable with 0.5. Clearly, w{B) decreases as B grows, signaling that 
the transition becomes stronger in the presence of the magnetic held. We will return to this observation 
below in Sec. 5. 

Besides being useful for quantifying the strength of the transition, the peak of the susceptibility 
allows for yet another determination of the transition temperature. Fitting for the peak maximum, 
we obtain Tc{x^^} = 113(4) MeV, consistent with the results obtained for all other observables. We 
mention that the peak positions for the up and down quark susceptibilities also agree within errors. 

Finally, in Fig. 5 we summarize our determinations of Tc in the QCD phase diagram. We consider 
the results for eB < 1 GeV^ obtained for the light quark condensates and for the strange quark number 
susceptibility [5]. In addition, we also include the transition temperatures at eB = 3.25 GeV^ obtained 
using the light quark condensates, the strange quark number susceptibility and the Polyakov loop. 
(Note that the inflection point of Pr at eB < 1 GeV^ is not pronounced enough to allow for a stable 
ht.) To interpolate Tc{'Pud} and ^c{c|} for all magnetic helds, we found the following function sufficient. 


Tc{eB) 


Tc{0) • 


1 + ai{eB)'^ 
l + a2(eH)2’ 


(3.1) 


giving the ht parameters shown in Tab. 1. The resulting ht is also shown in the hgure. 



eB (GeV2) 

Figure 5. The QCD phase diagram in the magnetic held-temperature plane. Previous results at weaker 
magnetic helds [5] are complemented by our hndings at high eB. The points have been slightly shifted horizontally 
for better visibility. The dotted and the dashed lines show an interpolation of the results for and for 

c|, respectively, according to Eq. (3.1). 



TM 

ai 

0-2 

^ud 

160(2) MeV 

0.54(2) 

0.82(2) 

^2 

174(2) MeV 

0.78(1) 

1.28(1) 


Table 1. Fit parameters of the function (3.1). 
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4 Results in the asymptotic magnetic field limit 


Our results in the right panel of Fig. 4 indicate that the transition becomes significantly sharper as 
the magnetic field increases. This observation raises the question: what happens if B is even larger? 
Does the crossover terminate and turn into a real phase transition? To answer this question, we have 
to consider the limit eB 3> Aqqq. Asymptotic freedom dictates that in this limit quarks and gluons 
decouple from each other. Still, the explicit breaking of rotational symmetry by B and the corresponding 
dimensional reduction in the quark sector [9] suggests that this limit is not simply given by a pure gluonic 
theory plus non-interacting (electrically charged) quarks. Indeed, based on the structure of the gluon 
propagator in strong magnetic fields. Ref. [25] has shown that the effective action describing this limit is 
an anisotropic pure gauge theory. The anisotropy amounts to an enhancement of the chromo-dielectric 
constant in the direction parallel to the magnetic field, characterized by the coefficient k, 


eB ^ Aqqd 




tr ;Bn -|- tr B\ -|- (1 -|- g^n{B)) tr £?< + tr 


k{B) oc B. (4.1) 


The definition of the gluonic field strength components B and £ is given in Eq. (A.2) below. The 
enhancement of the parallel chromo-dielectric constant implies that the corresponding field strength 
component Ty is suppressed. This tendency is already visible in our full QCD simulations at strong 
magnetic fields, see Fig. 11 in App. A below. 

Therefore, as B is increased, the QCD effective Lagrangian approaches the anisotropic gauge theory 
given by Eq. (4.1). Assuming that this theory has a first-order phase transition. Ref. [17] has conjectured 
that the strong magnetic field region of the QCD phase diagram should exhibit a critical point. Here we 
address this question in more detail. First of all, in App. A, we reproduce the results of Ref. [25] for the 
magnetic field-induced anisotropy using the effective action in the Schwinger proper-time formulation. 
The resulting anisotropic gauge theory can be simulated directly on the lattice. The setup and the 
simulation algorithm are described in App. B. The main difference to simple pure gauge theory amounts 
to multiplying the plaquettes lying in the z — t plane by the anisotropy coefficient k. The exact form 
of the anisotropic action is given in Eq. (B.2). 

Before discussing the lattice simulations of the anisotropic theory, let us make one more remark. 
Besides writing down the effective Lagrangian (4.1), Ref. [25] also predicted that the scale Aqcd of this 
theory (generated through dimensional transmutation) is much smaller than the QCD scale at H = 0: 
Aqcd Aqcd for a very broad range of magnetic fields. In the absence of further dimensionful scales in 
the anisotropic theory, this implies that the deconfinement transition temperature for strong mangetic 
fields is also much smaller than Tc{B = 0). Below we will also address this prediction. 

Due to the exact ^(3) symmetry of the anisotropic theory, the deconfinement transition is charac¬ 
terized by the projected Polyakov loop (2.6). In the left panel of Fig. 6, this observable is plotted against 
the inverse gauge coupling /3 for several values of k, as measured on the 16^ x 4 lattices. At k = 0, we 
reproduce the results of Refs. [30, 31] - in particular, the deconfinement transition occurs at /3c ~ 4.07. 
The results indicate f3c to be strongly reduced as k grows^. We find that f3c scales approximately with 
1 /\/k, see the right panel of Fig. 6. Extrapolating to k = oo we obtain /3 c(k = oo) ~ 2.42(5). Besides 
approaching this limit via finite values of the anisotropy coefficient, we also develop an algorithm to 
simulate directly at k = oo. The corresponding setup is described in App. B. 

'^Here we simulated at fixed values of the ratio k// 3. This continuous rescaling of k has no effect on, for example, the 
critical inverse coupling. 
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Figure 6. Left panel: the projected Polyakov loop as a function of the inverse gauge coupling for various values 
of the anisotropy coefficient, as measured on the 16^ x 4 lattices. The solid lines merely serve to guide the eye. 
Right panel: the critical inverse coupling as a function of 1/ y/n. The extrapolation to k = cx) is compared to 
the result of the direct simulation at infinite anisotropy (the latter shifted horizontally for better visibility). 


The left panel of Fig. 7 shows at infinite anisotropy. We find that the critical inverse coupling 
on the 16^ x 4 lattice is comparable with the extrapolation based on finite anisotropies, see the right 
panel of Fig. 6. In addition, a comparison of the results at different spatial volumes 12^ ... 24^ reveals 
that the transition becomes sharper as the volume increases, as typical for real phase transitions. We 
also repeated this analysis on Nt = 8 lattices, see the right panel of Fig. 7. The critical couplings 
are clearly different, showing that the transition is indeed related to the finite temperature. We also 
mention that finite volume effects in j3c are observed to be unusually large - above 10% for Nt = 4. 
(For comparison, the finite volume effects at k = 0 on similar lattices are of 0.1% [30].) We suspect 
that this is due to lattice artefacts - indeed, the effect is considerably smaller for Nt = 8, see the right 
panel of Fig. 7. 




Figure 7. The projected Polyakov loop as a function of the inverse gauge coupling at k = oo, for various 
lattice volumes with Nt = A (left panel) and Nt =8 (right panel). The solid lines merely serve to guide the eye. 


To determine the nature of the transition, we calculated the susceptibility of the projected Polyakov 
loop (2.7). This observable is shown in the left panel of Fig. 8, with a normalization by the spatial 
volume. Within statistical errors, the height of the normalized susceptibility peak is observed to be 
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2.4 2.5 2.6 2.7 



Figure 8. Left panel: the susceptibility of the projected Polyakov loop, normalized by the spatial volume, as 
a function of the inverse coupling. Various spatial volumes with TYj = 4 are compared. The solid lines serve to 
guide the eye. Right panel: histogram of near the critical temperature on the 16^ x 4 lattices. 


independent of V. In other words, the peak height scales linearly with V, which we take as strong 
evidence that the transition is of first order. The histogram of at fi = 2.4855 as measured on the 
16^ X 4 lattices is shown in the right panel of Fig. 8, revealing the two-peak structure characteristic for 
first-order transitions. 

Through the equivalence between this anisotropic gauge theory and QCD with asymptotically 
strong magnetic helds, the above hnding implies that the QCD phase diagram exhibits a critical point 
in the strong magnetic held region, where the crossover turns into a real phase transition. Based on our 
full QCD results for the light quark susceptibilities, we will estimate the magnetic held corresponding 
to the critical point in Sec. 5. 

The next step is to relate the critical inverse coupling fic to the critical temperature Tc in physical 
units. To do so, we must set the lattice scale /3(a). In principle, the magnetic held is not expected to 
change this scaling relation (cf. Ref. [5]). However, to arrive at our anisotropic gauge theory, B has been 
taken to inhnity, i.e. it also exceeds the squared lattice cutoff a“^. Clearly, the lattice scale determined 
at P = 0 becomes invalid beyond this point. Thus, in order to determine the lattice spacing, one needs 
a dimensionful quantity whose value is known in the asymptotic limit - for example a purely gluonic 
observable, where the P-dependence is expected to be only mild. A possible candidate for this role 
is the parameter wq dehned from the gradient how of the gauge links [32] that is often used for scale 
setting in QCD, as suggested in Ref. [33]. 

We determined wq on our zero-temperature full QCD ensembles [6] for eP < 1 GeV^ and also for 
eP = 3.25 GeV^ at our lowest temperature T ~ 75 MeV. The results for the ratio wo/wo{B = 0) are 
plotted in the left panel of Fig. 9, showing a mild reduction of this parameter as P grows. A ht of the 
form similar to Eq. (3.1) describes the data well and suggests a saturation towards the asymptotically 
strong magnetic held limit. Nevertheless, we cannot exclude a signihcant dependence of wq on P for 
P > 3.25 GeV^. 

In addition, we can also gain some insight by considering the dimensionless combination TcWq. How 
close full QCD at eP = 3.25 GeV^ is to the asymptotic limit can then be quantihed by matching T^wq 
with the anisotropic theory. Multiplying our full QCD results for wq by the transition temperature (here 
we take the dehnition of Tc employing the inhection point of the strange quark number susceptibility, 
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cf. Fig. 5), TcWq is shown in the right panel of Fig. 9. Motivated by the scaling of /3c (cf. the right panel 
of Fig. 6 ), the results are plotted against Employing the result for wq from Ref. [33], at zero 

magnetic field we have Tc{B = 0) • wo{B = 0) = 0.174(3) GeV • 0.1755(19) fm = 0.155(3). 




eB (GeV2) (eB) (GeV ') 

Figure 9. Left panel: magnetic field-dependence of the parameter wq using various lattice spacings and a fit 
(dashed line) of the form similar to Eq. (3.1). Right panel: the dimensionless combination TcWq in full QCD 
(\l\feB > 0.5 GeV“^) and in the anisotropic pure gauge theory {IjyfeB = 0). The dashed line indicates the 
B = 0 limit. 

To carry out the comparison to the asymptotic limit, we also determined Wo/a on symmetric 16^ and 
24^ anisotropic gauge configurations at the critical couplings corresponding to the 16^ x 4 (Pc ~ 2.47) 
and to the 24^ x 8 (/3c ~ 2.7) lattices^. We observe that the combination TcWq = wo{f3c)/a ■ 1/Nt 
- similarly to /3c - suffers from large lattice discretization effects and exhibits a downwards trend 
towards the continuum limit. We take the result for the Nt = 8 data as an upper limit, giving 
limB^oo{TcWo) < 0.076. This value is also included in the right panel of Fig. 9. Altogether, the 
results are compatible with a monotonous dependence of TcWq on B. To extrapolate TcWq reliably to 
the continuum limit in the anisotropic theory requires further simulations on finer lattices and will be 
discussed in a forthcoming publication. 

To summarize, the lattice results favor a saturation of wq and a monotonous reduction of TcWq as 
the limit R —)> oo is approached. This suggests a monotonous reduction of Tc{B) towards the asymptotic 
limit. Nevertheless, based on the available findings, no final statement about lim 5 _>.oo 7c can be made. 

Let us make one more remark about the k = oo anisotropic theory. Since the parallel chromoelectric 
component of the action vanishes, all plaquettes lying in the z — t plane are unity. This implies 
that all Wilson loops W in this plane are trivial, and the static quark-antiquark potential oc logLF is 
independent of the distance. Accordingly, there is no force acting on quark-antiquark pairs if they are 
separated in the direction of the magnetic field, i.e. the string tension (T|| in this direction vanishes. This 
is in line with recent lattice determinations of the string tension in magnetic fields [34]. 

5 Conclusions 

In this paper, we determined the nature and the characteristic temperature of the chiral/deconfinement 
transition of QCD at an extremely strong background magnetic field eB = 3.25 GeV^. The results for 

^ Just as in full QCD, the gauge links are evolved here using the symmetric gradient flow. 
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various observables consistently show that the transition temperature is further decreased compared to 
its value at lower magnetic fields. For the light quark condensates, the reduction of Tc is due to the 
so-called inverse magnetic catalysis: between eB = 1 GeV^ and eB = 3.25 GeV^, and of are 
significantly reduced in the transition region. At the same time, the condensates are enhanced by B 
both for T <^Tc and for T ^ Tc (the latter effect is small, since the condensate is suppressed at high 
temperatures). 

Gomparing the behavior of the up and down quark condensates and that of the Polyakov loop also 
revealed that there is no splitting between the transition temperatures for the individual flavors, neither 
is there significant difference between the chiral and the deconfinement transition temperatures. On the 
contrary, the different definitions of Tc tend to approach each other as B grows and at eB = 3.25 GeV^ 
all observables exhibit a single transition temperature of around 109—112 MeV, see Fig. 5. Furthermore, 
we performed a finite size scaling analysis of the light quark susceptibilities, which has revealed that 
there is no singularity in the infinite volume limit and, thus, the transition remains an analytic crossover 
even at eB = 3.25 GeV^. 

In addition, we considered the asymptotically strong magnetic field limit, and simulated the cor¬ 
responding effective theory on the lattice. This limiting effective theory - an anisotropic pure gauge 
theory - was found to exhibit a first-order deconfinement phase transition. Together with our findings 
above, this implies the existence of a critical point in the QGD phase diagram. To provide a first 
estimate for the magnetic field i?cP corresponding to the critical point, let us return to our results 
about the width 'w{B) of the light quark susceptibilities. We have seen that the width is reduced as 
the magnetic field grows, see the right panel of Fig. 4. Assuming a linear dependence of w on B and 
extrapolating in the magnetic field we find that w vanishes at 

eBcP ~ 10(2) GeV^. (5.1) 

In light of the fact that the ^-dependence of some of our observables (e.g. of Tc and of wq) tends to 
flatten out as B grows, this first estimate should rather be taken as a lower bound for Bqp- We mention 
that in order to simulate with magnetic fields of strengths comparable to that in Eq. (5.1), lattices with 
Nt > 28 are required, out of reach for current computational resources. 

In the absence of a priori known dimensionful scales in the B ^ oo system, we could not determine 
lim^^oo Tc in physical units. Nevertheless, the deconfinement transition temperature of the anisotropic 
theory is expected to be much smaller than Tc{B = 0) [25]^. Our results for the combination TcWq 
are compatible with this prediction and suggest a gradual reduction of the deconfinement transition 
temperature as B is increased^. Taking these aspects into account. Fig. 10 represents a sketch of the 
deconfinement transition line in the QGD phase diagram for a broad range of magnetic fields. 

The reader might wonder whether it is possible that the crossover at eB < 3.25 GeV^ and the 
first-order transition in the asymptotic limit are not connected by a single line. To see that this is not 

^ The discussion in Ref. [25] bases on renormalization group arguments and on the separation of scales Aqcd ^ rn^ <C 
VeB (here Aqcd is the dynamical scale of the large-B theory) to conclude that Aqcd <C Aqcd for a very broad range 
of magnetic fields - in fact, up to eB being millions of orders of magnitudes larger than Aqqq. To find the complete 
B-dependence of the running of the strong coupling and to prove rigorously that Aqcd ^ Aqcd even in the B —>■ oo 
limit, a full treatment of the divergent one-loop Feynman diagrams in the presence of background magnetic fields would 
be necessary. Without relying on the lowest-Landau-level approximation - which might not be justified for the case of 
divergent diagrams - this is a very difficult task. 

Note that our setup at B —>■ oo describes the low-energy effective action for gluons. Thus, the results in the anisotropic 
theory have no implications for the chiral transition. For more details on this point, see Ref. [25]. 
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Figure 10. The deconfinement transition temperature against the background magnetic field. The results 
of our full lattice QCD simulations (white background) are complemented by the prediction (gray background) 
based on the results corresponding to the B ^ oo limit and on the extrapolation of the light quark susceptibility 
peak to high magnetic fields (see the text). 

the case, note that by varying the anisotropy parameter k, one can continuously deform the anisotropic 
theory to usual pure gauge theory, as was demonstrated in Fig. 6. Furthermore, the isotropic pure 
gauge theory can be thought of as QCD with infinitely heavy quarks and thus can be continuously 
transformed into full QCD by increasing the inverse quark masses from zero to their physical values. 
Thus, the transition we identified at i? —)> oo is indeed the same deconfinement transition that occurs 
at low magnetic fields. 

Let us highlight that according to this discussion, having a decreasing deconfinement transition 
temperature is actually natural to QCD. Furthermore, since the B ^ oo limit is independent of the 
quark masses®, a similar reduction of Tc by the magnetic field should also take place in QCD with 
heavier-than-physical quarks. However, in the latter case this reduction most probably follows an 
initial increase in the transition temperature, cf. Refs. [5, 35]. Indeed, recent lattice results employing 
overlap fermions and pion masses of about 500 MeV indicate inverse catalysis to occur around the 
transition temperature at the magnetic field eB ~ 1.3 GeV^ [8]. 

Finally, we note that magnetic fields well above the strength (5.1) are predicted to be generated 
during the electroweak phase transition in the early universe [36]. If these fields remain strong enough 
until the QCD epoch, the emerging first-order phase transition might have several exciting consequences. 
Via supercooling, bubbles of the confined phase can be formed as the temperature drops below T^, 
leading to large inhomogeneities, important for nucleosynthesis [37]. Collisions between the bubbles 
can also lead to the emission of gravitational waves and, thus, leave an imprint on the primordial 
gravitational spectrum [38]. An absence of such signals, in turn, would imply an upper limit for the 
strength of the primordial magnetic fields. 
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A Effective action in the asymptotic magnetic field limit 


in this appendix we demonstrate how asymptotically strong magnetic helds induce an anisotropy in the 
gluonic sector using an Euler-Heisenberg-type approach. The QCD effective Lagrangian in Euclidean 
space-time is 

= + C^B,G^,) = - logdet + m/] , (A.l) 

f=u^d^s 

and the quark determinant will be regularized using Schwinger’s proper time formulation [39]. Since 
the electromagnetic held exceeds all scales in the system and in particular, (eB)'^ ^ trG^^, we may 
approximate the the chromo-helds in the fermionic action to be weak. In addition, we assume the 
chromo-helds to be covariantly constant, D^G^p = 0 to enable a fully analytical treatment of the 
problem. Given this condition, the held strength can be gauge transformed to be constant in space- 
time and diagonal in color space [40], Gpy = d\&g{Gpyc) with the color index c = 1, 2, 3. 

Let us decompose the chromo-helds to chromomagnetic/chromoelectric components, 

= = + (A.2) 

parallel or perpendicular to the electromagnetic held B. The leading terms in the strong i?-expansion 
are quadratic in the chromo-helds and thus, to hnd the coefficients of the respective components, it 
suffices to consider separately the effect of B and i3||, B and B±, B and Ty and B and The effective 
Lagrangian for these components for small background magnetic helds was determined in Refs. [29, 41]. 
A similar calculation, generalized to hnite temperatures and constant Polyakov loop backgrounds was 
performed in Refs. [7, 42]. 

Let us hrst take the case of B and i3||. For each havor we may choose our coordinate system such 
that qfB is positive. Then, each color component experiences a total (positive) magnetic held qfB + B^^c 
so that 

1 V—A n /" ds {QfB B\\„)s 

+ ^coth - -2 -. (A.3) 

f,c f 

Since B^^^ only appears in the sum with q/B, the effective Lagrangian becomes independent of the 
chromomagnetic held in the limit qjB 3> B\\f.- This implies that quarks become insensitive to Sy, i.e. 
decouple from this gluonic component. 

Next we take the case with B and B_i. Rotating our coordinate axes for each color component such 
that the z axis points in the direction of the total magnetic held we get 


C‘^{B,B^) = 


1 


E 

f,C 


m. 


{q,Bf + 


ds , sjivB? + Si,s 

coth - - 


ml 


(A.4) 


In the strong B limit, this becomes independent of B±c, signaling that quarks decouple from the 
perpendicular chromomagnetic component of the gluons as well. 
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For a perpendicular chromoelectric field, for each color component we can perform the Lorentz 
transformation that eliminates the electric field and, in turn, gives a total magnetic field ^ 

The corresponding effective Lagrangian equals Eq. (A.4), but with B±c replaced by £±c- This implies 
the decoupling of quarks from the perpendicular chromoelectric fields. 

Finally, for a parallel chromoelectric field Ty, we have a Landau problem in the x — y as well as in 
the z — t planes, giving 





(A.5) 


Taking the limit qfB ^ T||(,,mp we see that - unlike for the other components above - a non-trivial 
dependence on Ty remains. We are interested in the quadratic term, proportional to trT|, which 
contributes® to the gluonic Lagrangian tr of Eq. (A.l), 


C'i{B,0{8l)) 



k{B)^ 


1 

24 ^ 


k//e 

/ 


m) • 


(A.6) 


Altogether, the asymptotically strong magnetic field limit of the QCD effective Lagrangian indeed 
equals Eq. (4.1). Thus we find that the chromo-dielectric constant is enhanced in the direction of the 
background magnetic field, and the coefficient k,{B) coincides with the result of Ref. [25]. 

Having k ;§> 1 in the action implies that 
the corresponding gluonic field strength com¬ 
ponent trT| is strongly suppressed. In other 
words, the anisotropy in the chromoelectric 
part of the action density, 

^(^) = -trTf^, (A.7) 

is enhanced by H. To back up this prediction, 
in Fig. 11 we plot A{£) as a function of the 
magnetic field, based on our zero-temperature 
results at eB < 1 GeV^ [29] and the measure¬ 
ments at eB = 3.25 GeV^ at our lowest tem¬ 
perature T ~ 75 MeV. The results clearly 
indicate that the anisotropy is positive and 
strongly increased as B grows. We note that 
due to the coupling between the gluonic field strength components, a similar anisotropy in the chromo- 
magnetic sector also appears, altogether giving rise to the hierarchy trHy > trH^ = trT^ > trT| at 

® Here we omitted a divergent term of the form trf|^ log A/m/, where A is a cutoff entering as the lower endpoint of 
the proper time integration sq oc 1/A^. This divergence can be eliminated by the multiplicative renormalization of the 
wave function f|| and of the gauge coupling g [39]. Closer inspection of Eqs. (A.3) and (A.4) shows that the same type of 
divergence is present for the other components as well. Thus, these B-independent terms merely represent an isotropic 
redefinition of the gauge coupling g, which does not alter the form of the effective Lagrangian for strong magnetic fields. 
Another divergence, independent of the gluonic field strengths, takes the form (q/B)^ log A/m/ and is canceled by the 
renormalization of B and of q/ [39]. Thus, the necessary renormalizations at B —>■ (xj are of the same type as for the 
theory at small magnetic helds. 



Figure 11. The anisotropy of the chromoelectric gluonic 
field strength component in full QCD. 
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low temperatures. The same hierarchy is also observed in the anisotropic gauge theory^. 

We note that the calculation leading to Eq. (4.1) can also be performed for nonzero temperatures. 
At T > 0 an additional factor appears in the proper time integral due to the sum over Matsubara 
frequencies, containing an elliptic ©-function. This factor decouples from the S-dependence, implying 
that even for T > 0, only the chromo-dielectric constant is affected. The coefficient k is, however, 
altered as 


k{B,T) 


1 

247r2 


f 


m} 


dse ©3 


g-m2/(4sT2) 


(A.8) 


The integral over s equals unity at T = 0 and is reduced monotonously (and smoothly) as the tempera¬ 
ture grows. Simulating the anisotropic gauge theory according to the Lagrangian (4.1) on the lattice, we 
found that the theory exhibits a first-order phase transition. Thus, since the smooth k{T) dependence 
does not affect the discontinuous transition, in order to locate the critical temperature it suffices to 
simulate the theory at fixed (large) k values. 


B Simulating anisotropic pure gauge theory on the lattice 

In this appendix we discuss the simulation algorithm for the anisotropic pure gauge theory described 
by the Lagrangian (4.1). The corresponding path integral 

Z = I VU (B.l) 

can be simulated directly on the lattice. Here, U denotes the gauge links, /3 = Q /is the inverse gauge 
coupling and the anisotropic gauge action reads 

^amso ^ ^ 1 Pg I-J. 

IJ,<U 

where are linear combinations of closed loops lying in the g — v plane. We take the tree-level 
Symanzik improved gauge action such that these loops include the 1x1 plaquettes ^ and the 2x1 
rectangles ^ with appropriately tuned coefficients [43], 

P,u = -^ (1 - Ul"") + 3 (1 - ■ (B.3) 

The correspondence between the continuum and lattice expressions reads 

tr = 2 Re tr Ptz , tr = Re tr [Ptx + Pty ], tr = 2 Re tr P^y , tr Bj_ = Re tr [Pyz -t- Pxz] ■ 

(B.4) 

To simulate this theory, we use an overrelaxation/heatbath algorithm, based on the isotropic pure 
gauge implementation by the MILC collaboration [44]. One trajectory consists of one overrelaxation 
step followed by four heatbath steps. The simulation at finite anisotropy coefficient k simply involves 
multiplying the plaquettes and rectangles lying in the z — t plane by k. We observe that autocorrelation 

^ To see how the anisotropic dielectric constant affects the chromomagnetic components, it is instructive to consider the 
gauge potential A large value of k implies a suppression of tr£’|^ and a corresponding suppression of the fluctuations 
in Az and in At. This suppression propagates into the magnetic sector and creates the anisotropy between trS^. ^rid 
irB?!. Indeed, while the former contains Az, the latter does not. 


1 + k{B)//3, IJ, = Z,l' = t, 

1, otherwise, 


(B.2) 
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times grow large as k increases, similarly to the issue of critical slowing down of the isotropic theory 
at large /?. This prohibits approaching k —)• oo, necessary for the asymptotically strong magnetic field 
limit. However, it is possible to modify the algorithm to simulate directly at k = oo. In this limit, the 
z — t component of the action and the remaining five components decouple, and the links are restricted 
to the subspace H[?7] of configurations, where P^t is minimal. This subspace is defined by 

n[U] = {U^\Ul,^^ = l}. (B.5) 

Indeed, any fluctuation in the link variables that leads off of this subspace makes the action infinitely 
large and is thus forbidden. (Note that if equals the unit matrix, then so does 

We thus have to parameterize the subspace H[?7] in terms of the gauge links I/^. Let us label the 
lattice sites by n = (na,, n^, n^, n*) with 0 < < N^. To find the parameterization of Eq. (B.5), it is 

advantageous to fix the links to 1 on a so-called maximal tree. The specific choice for the tree is shown 
in the left panel of Fig. 12. (Note that Faddeev-Popov fields are absent for such a gauge fixing [45].) 
In order to have unit plaquettes for Hz < Nz — 1 and m < Nt — 1, all t-links must be set to unity at 
these sites. To have unit plaquettes on the last z-slice, all z-links at Nz — 1 must be set equal, denoted 
by Lz- Similarly, all the t-links at Nt — 1 must be set equal, denoted by Lf, see the visualization in the 
right panel of Fig. 12. These remaining links correspond to the local Polyakov loops in the z- and in 
the t-direction, lying in the z — t plane at a given and Uy. Finally, to ensure that the plaquette at 
the corner Uz = Nz — 1, nt = Nt — 1 is unity, we need LzLtL\L\ = 1, i.e. Lz and Lt must commute. 
Altogether, the subspace in question reads 

L![[/] = {U^ I Uz{n) = \ynz^Nz- 1, 

Ut{n) = 1 Vn* / W - 1 , 

Uzipx^ nyi Nz 1 , nt) — Lzijixt ^y) (^^' 6 ) 

Ut{nx,ny, nz,Nt - 1 ) = Ltirix, ny) tfnz, 

[Lzijix^ ny), Lti^nxj ny)^ — 0 }. 

Notice that the ‘degenerate’ timelike Polyakov loop Lt (represented by the blue arrows in the right panel 
of Fig. 12) appears multiple times in the action - in fact, in ANz plaquettes and in 12Nz rectangles. 
The corresponding ‘staples’ are all taken into account in the update of Lt (and similarly for Lz). 

There are several ways to fulfill the commutativity relation [Lz,Lt\ = 0. One possibility (setup A) 
is to simply set = 1 for all and ny. Another approach (setup B) is to constrain Lz to be a center 
element, Lz{nx,ny) G Z 3 . The two setups only differ on a set whose measure vanishes in the limit, 
where all lattice extents are taken to infinity. Note that the expectation value of the average z-Polyakov 
loop [defined similarly as the usual Polyakov loop P, Eq. (2.5)] is three for setup A, whereas it 
is zero for setup B, if the spatial size of the system is large enough. Nevertheless, we checked that 
observables sensitive to the finite temperature transition (P, the gauge action, etc.) all have vanishing 
correlators with P^^'i . In fact, we found that the setups A and B give identical results for P and for 
the gauge action for all values of the inverse gauge coupling f3 on the 16^ x 4 lattices. In other words, 
center symmetry breaking in the z direction appears to be completely irrelevant for the deconfinement 
phase transition. 

In addition, we also tried allowing both Lz and Lt to be general SU(3) matrices (which violates 
the commutativity relation). This approach (setup C) turned out to introduce negligible differences in 
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Figure 12. One z — t plane of the lattice with periodic boundary conditions in both directions. The dashed 
arrows at Uz = Nz {jit = Nt) indicate the copies of the links at riz = 0 {rit = 0). Left panel: the gauge links 
on the maximal tree (yellow arrows) are fixed to the unit matrix. Right panel: in order for the configuration to 
belong to the subspace (B.6), further gauge links are set to unity (yellow) and the green and blue links are set 
equal, respectively. Finally, the blue and green ‘degenerate’ Polyakov loops must commute with each other. 


the results®. For the simulations presented in the body of the paper, we considered setup A and set 
Lz = 1 throughout the lattice. As a consistency check, besides the overrelaxation/heatbath algorithm, 
we also considered a hybrid Monte-Carlo update and found that the two give fully consistent results. 
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